Comment on the use of the method of images for calculating electromagnetic 

responses of interacting spheres 
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In this Comment I argue that the method of images used by Huang, Yu, Gu and co-authors [Phys 
Rev. E, 65, 21401 (2002); Phys Rev. E., 69, 51402 (2004)] to calculate electromagnetic properties of 
interacting spheres at finite frequencies is inapplicable and does not provide a physically meaningful 
approximation. 



Recently, Huang, Yu, Gu and co-authors (referred to 
as HYG below) have applied the method of images to 
study theoretically the electromagnetic properties of two 
interacting spherical particles 0,0- As is well known, 
the method of images can be applied to spherical con- 
ductors in the electrostatic limit, i.e., when the dielectric 
constant e can be formally set to ioo and the Bergman- 
Milton spectral parameter s = l/(e — 1) is equal to zero. 
At finite frequencies, when s is not small compared to the 
generalized depolarization factors s n , the method of im- 
ages is not applicable. However, HYG apply the method 
to dielectric particles at arbitrary frequencies, assuming 
only that the size of the two-sphere dimer is much smaller 
than the external wavelength. In particular, they claim 
to be able to extract the factors s n and the correspond- 
ing oscillator strengths F n , which characterize the elec- 
tromagnetic response of a system within the quasistatics. 
In the first paper of the series |3( and in Rcf. jlj the au- 
thors mention that their method is approximate. How- 
ever, in the more recent paper Q it is presented as exact 
and used without restriction. In the present Comment 
I show that it is impossible to calculate the quantities 
s n and F n using the method of images. Moreover, the 
expressions for s n and F n given by HYG are not con- 
sistent with the exact electrostatic solution. Thus, the 
mathematical formalism developed by HYG is not only 
not exact, but does not provide a physically meaningful 
and controllable approximation. 

We start with a brief review of mathematical formal- 
ism used by HYG. Within the quasistatics, dipole mo- 
ment of an arbitrary particle characterized by the dielec- 
tric function e(u) and excited by a homogeneous external 
field Eq exp(— iut) can be written as <iexp(— itut) where 
d = 6iE{). Here a is the polarizability tensor. If polariza- 
tion of the external field coincides with one of the prin- 
cipal axes of a, both vectors d and Eq become collinear. 
The corresponding scalar polarizability can be written in 
the Bergman-Milton spectral representation £| as 



polarizability tensor coincides with the axis of symmetry, 
and the other two axes are perpendicular to the first one 
and to each other, but otherwise arbitrary. HYG consider 
two interacting spheres of the radius R each separated by 
the center-to-center distance 2L, obtain the diagonal ele- 
ments of the polarizability tensor and derive the following 
expressions for F n and s n : 



F (L) _ F (T) _ p _ 

4n(n+ l)sinh 3 aexp[-(2n + l)a] ,(2) 
#> = i{l-2exp[-(2n+l)a]} , (3) 



>i T) = i{l + exp[-(2n+l)a]} , 



(4) 



n= 1,2,3,... 



where the upper index (L) denotes longitudinal modes, 
(T) denotes transverse modes, and a is the solution 
to cosha = L/R, or, explicitly, a = \n[L/R + 
J (L/R) 2 — 1] It can be verified that F n satisfy the 
sum rule Y,n F n = l- 

Everywhere below we consider only the longitudinal 
modes, although the results of HYG for the transverse 
modes are also incorrect. The longitudinal modes are 
more important physically, since they are known to 
produce extremely high field enhancements in axially- 
symmetrical arrays of nanospheres |fij and have been ex- 
tensively studied in conjunction with the single-molecule 
spectroscopy Q. 

First, let us discuss the small- frequency limit for con- 
ductors. In this limit, the dielectric function can be writ- 
ten as e = Aitia/uj, where a is the static conductivity. 
Correspondingly, s cx uj — > 0, and we can expand a into 
a power series in s. The expansion can be obtained from 
||TJ and reads 



4vr ^ 



(1) 



where v is the volume of the particle, s n - the generalized 
depolarization factors satisfying < s n < 1 and F n are 
the corresponding oscillator strengths. 

In the case of two spheres, one principal axis of the 



fc=0 

A k =J2F n /4 +1 



(5) 
(6) 



The electrostatic polarizability is given by a cs = 
(v/4tt)Aq. The method of images can provide an exact 
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expression for a cs and, correspondingly, for Aq. How- 
ever, since there is an infinite number of different sets of 
F n , s n that produce the same value of Aq, it is impossible 
to find these coefficients from the electrostatic solution. 
We emphasize that this is not possible even if one con- 
siders the inter-sphere separation as an additional degree 
of freedom, since all quantities (Ak,F n and s n ) depend 
parametrically on L/R. Instead, if the summation in the 
right-hand side of © is truncated at n = N, one needs 
to calculate all coefficients Ak from k = to k — 2N — 1 
in order to make the system of equations © sufficiently 
determined. But the electrostatic solution based on the 
method of images can provide only one of these coeffi- 
cients, namely, Aq. 
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Fig. 1. Ao(L/R) — ^4o(oo) as a function of the relative 
separation L/R calculated by different methods. 

Although one can not expect that the set of F n ,s n 
given by HYG ©7© would produce, upon substitu- 
tion into 10, the correct expansion coefficients Ak for 
k > 0, it is still possible that the value of Aq ob- 
tained in this manner is correct. However, as is demon- 
strated in Fig. 1, this is not so. In this figure, we plot 
the function A (L/R) — A (oo) (for conducting spheres, 
Aq(oo) = 3) calculated by different methods. The math- 
ematically rigorous result is shown by the solid curve, 
and the result of HYG by the long dash. We also show 
in this figure two analytical asymptotes valid for L » R 
(shorter dash) and L — > R (dots). The different curves 
in Fig. 1 are explained below in more detail. At this 
point, we note that the result of HYG for Aq(L/R) is ac- 
curate at large separations (L ^> R), but breaks down 
when L/R ps 1.2, and becomes grossly inaccurate at 
L/R 1.03. In particular, the HYG curve has a sin- 
gularity at L/R = x c ee (2 2 / 3 + l)/2 4 / 3 w 1.026. This 
is due to the fact that the first depolarization factor s\ 
defined by © crosses zero when L / R = x c . The appear- 
ance of negative depolarization factors for smaller inter- 
sphere separations is unphysical and can, in particular, 
result in divergence of the electrostatic polarizability |8| . 

In the next two paragraphs I explain how the data 
for different curves shown in Fig. 1 were calculated. The 
solid curve was obtained by diagonalization of the electro- 
magnetic interaction operator W whose matrix elements 
are given [jj by 
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(i + iy. 



l+l' 



IV 



(21 + l)(2l' + 1) (L/RY+^+HU'l ' 



(7) 



where i, i' — 1,2 label the spheres, 1,1' = 1,2,... and Zi 
is the z coordinate of the center of ith sphere, assuming 
the z-axis coincides with the axis of symmetry. The de- 
polarization factors s n are the eigenvalues of W while the 
oscillator strengths can be found as squared projections 
of the corresponding eigenvectors \n) on the vector of ex- 
ternal field: F n = (E\n)(n\E), where \E) is normalized 
so that (E\E) = 1 8J. The matrix defined in Q) was 
truncated so that Z, V < Z max = 1000 and diagonalized 
numerically. In the absence of round-off errors and in 
the limit Z max — > 00, such diagonalization would produce 
the infinite set of exact values s„, F n . We note that at 
l max = 1000 and L/R> 1.01, all the modes whose oscil- 
lator strength are not very small (i.e., greater than 0.001) 
have converged with a very high precision, and that the 
round-off errors do not influence the results in any no- 
ticeable way since the matrix W is well-conditioned. 

The dots and short dash in Fig. 1 show the theo- 
retical asymptotes obtained by Mazets who has derived 
an ex pres sion for Aq in terms of hypergeometrical func- 
tions [lOj. He has also provided simple asymptotic ex- 
pansions which are valid for small and large inter-sphere 
separations. Thus, for longitudinal excitations, 
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(9) 



where Q(x) is the Riemann zeta- function and C is the 
Euler constant. The second term in the right-hand side 
of 10 is a correction due to the dipole-dipole interaction 
while the third term describe the next non-vanishing in- 
put due to the higher multipole interaction. It can be 
verified that the asymptotic expansion of the HYG re- 
sult coincides with © at least up to the sixth order in 
L/R. However, the small-separation asymptote iJSJ is 
dramatically different from the one that follows from the 
HYG formulas. 



3 



0.8 
0.6 - 
0.4 - 
0.2 - 





-q— r 



Exact -■&■ 
HYG --B- 



B- 



L/R = 1.2 



■0, \ 



_|_ 



0.2 



0.8 







0.6 - 
0.4 - 
0.2 - Q- 



0.3 0.4 
1 



0.5 



0.6 



-B... 



\ 



L/i? = 1.1 



-I — 1 — h'^ai ^e- 



0.1 



0.2 



F n 

0.4 - 



0.2 - 



i 1 1 r 

o- . _ 

~-<5 



■B- 



0.3 0.4 0.5 

~i r 



0.6 



J L 



L/R= 1.05 



I n ii I" i h i i ' i ^ I 



-e- 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 



0.3 



0.2 



0.1 - 



-Ox 



1 

L/R = 1.01 



■ B 



B'B 



...-B" 




I Ml i 



-0.2 







0.2 



0.4 



0.6 



Fig. 2 Bergman-Milton depolarization factors, s„, and the 
corresponding oscillator strengths, F n , for different relative 
inter-sphere separations L/R. Dashed lines are plotted to 
guide the eye. 
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Fig. 3. Dimensionless extinction parameter e = a e /kv as a 
function of wavelength A, where a e is the extinction cross sec- 
tion, k — 2ir/\, v is the total volume of the scatterer, plotted 
for different relative inter-sphere separations L/R. Polar- 
ization of the incident field is parallel to the axis of symmetry. 



Next, we compare the coefficients s n ,F n defined by 
(0) , J3J according to HYG with respective values obtained 
by direct diagonalization of the interaction matrix W. 
The results are shown in Fig. 2. A significant discrepancy 
already exists at L/R — 1.2 and becomes more dramatic 
as this ratio approaches unity. Negative depolarization 
factors are present in the plot for L/R = 1.01. We note 
that the smallest inter-sphere separation considered by 
HYG was L/R=l + 1/30 w 1.033. As was mentioned 



above, the negative depolarization factors appear for 
L/R < x c w 1.026. At these separations, results of any 
calculation based on the HYG formalism are expected to 
be grossly inaccurate and unphysical. However, this fact 
is not explained in Refs. For example, the choice 

of values for L/R in Fig. 5 of Ref. appears to be ran- 
dom, while, in fact, all these values satisfy the critical 
condition L/R > x c . 

While it is demonstrated in Fig. 2 that the values of 
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s n , F n calculated according to HYG are inaccurate, these 
coefficients are not directly measurable in an experiment. 
However, they can be used to calculate various physically 
measurable quantities. For example, the extinction cross 
section is given by a e = Airkvlm J2 n Fn/ (s+s n ). In Fig. 3 
we plot the extinction spectra of two silver nanoparticles 
obtained for the same inter-sphere separations as in Fig. 2 
and for the longitudinal polarization of the external field. 
Interpolated data for silver from Ref. have been used 
to calculate the spectral parameter s as a function of 
wavelength. It can be seen that the spectra calculated 
using the formulas f° r s n ,F n differ dramatically 

from those calculated with the use of exact values of these 
coefficients. The discrepancy is evident even at relatively 
large separation, L/R = 1.2. It should be noted that in 
the case L/R= 1.01 the HYG spectra exhibit unlimited 
growth with the wavelength which starts in the near-IR 
region (data not shown). This is due to the appearance of 
negative depolarization factors and contradicts the gen- 
eral sum rules for extinction spectra which imply that a e 
must decrease faster than 1/A in the limit A — > oo @. 
Note that the presence of negative depolarization factors 
can result in even more severe anomalies of extinction 
spectra in dielectrics whose static dielectric permeability 
is positive, as well as the value of s in the limit A — ► oo. 

The papers 0, contain a number of other less sig- 
nificant inaccuracies. In particular, HYG confuse ori- 
entational averaging (for randomly-oriented bispheres) 
with the averaging over polarization of the incident light. 
Thus, for example, Eq. 2 in Ref. [l| is presented as a re- 
sult of averaging over polarization for a fixed bisphere. 



However, such averaging should clearly depend on the di- 
rection of the incident wave vector relatively to the axis 
of symmetry of the bisphere. In fact, the first equal- 
ity in this formula gives the result of orientational av- 
eraging, except that HYG are mistaken in stating that 
(cos 2 9) = (sin 0) = 1/2. It is easy to check that 
(cos 2 9) = 1/3 and (sin 2 9) = 2/3. Note that the second 
equality in Eq. 2 of Ref. [l| would be correct if the aver- 
aging is done over polarizations of the incident beam for 
a fixed bisphere, assuming that the incident wave vector 
is perpendicular to the axis of symmetry. 

It should be noted that on page 4 of Ref. 0, the au- 
thors acknowledge that the method of images is only ap- 
proximate but state that the approximation is very good 
and make a reference to the earlier work Q to support 
that statement. However, in Ref. [3j verification of the 
accuracy of the method of images is only done for rela- 
tively large separations, namely L/R > 1.5 (Figs. 3, 4 
in Ref. At these separations, the multipole effects 

are generally not important, which clearly follows from 
the data shown in these figures. However, in later pub- 
lications, HYG have used the method for much smaller 
separations, typically, L/R= 1 + 1/30. 

Finally, on the same page of Ref. [j], the authors write: 
"More accurate calculations based on bispherical coordi- 
nates can be attempted." This was, in fact, done in the 
above-referenced paper by Mazets [lCj , although only for 
perfect conductors. More general analytical results can 
be obtained with the use of the theory of hypercom- 
plex variables (a generalization of the conformal map- 
ping) [12. 
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